package com.gitee.wsl.mathematics.geometry

import com.gitee.wsl.ext.number.approx
import com.gitee.wsl.mathematics.function.JenkinsTraub
import com.gitee.wsl.mathematics.geometry.Polynomial.Companion.One
import com.gitee.wsl.mathematics.geometry.Polynomial.Companion.Zero
import com.gitee.wsl.mathematics.statistics.Functions.binomial
import kotlin.math.abs
import kotlin.math.absoluteValue
import kotlin.math.max
import kotlin.math.min
import kotlin.math.pow
import kotlin.math.sign
import kotlin.math.sqrt

class NoMoreDerivativeException(msg: String) : Exception(msg)

class Polynomial(val coefficients: List<Double>) {
    companion object {
        fun fromRoots(vararg roots: Double) =
            roots.map { Polynomial(1.0, -it) }.fold(One) { v, cur -> v * cur }

        val Zero = Polynomial(0.0)

        val One = Polynomial(1.0)

        val X = Polynomial(1.0, 0.0)

        val OneMinusX = Polynomial(-1.0, 1.0)

        fun Num(n: Double) = Polynomial(n)
    }

    constructor(vararg coefficients: Double) : this(coefficients.asList())

    inline val n get() = coefficients.size - 1

    inline val isEmpty get() = n == -1

    inline val isNumber get() = n == 0

    inline val isLine get() = n == 1

    inline val isSimple get() = coefficients.firstOrNull()?.approx(0.0) == false

    fun eval(x: Double) = coefficients.fold(0.0) { r, c -> r * x + c }

    fun eval(x: Polynomial) = coefficients.fold(Zero) { r, c -> r * x + c }

    val derivative by lazy {
        if (n < 0) {
            throw NoMoreDerivativeException("polynomial has no derivative")
        } else {
            Polynomial(List(n) { coefficients[it] * (n - it) })
        }
    }

    fun derivativeOrNull() = if (n < 1) null else derivative

    fun simplify(): Polynomial {
        val a = coefficients.find { !(it approx 0.0) } ?: return Zero
        val simplified = Polynomial(coefficients.dropWhile { it approx 0.0 }.map { it / a })
        return if (simplified.isEmpty) Zero else simplified
    }

    val roots by lazy { solve(coefficients) }

    override fun toString(): String {
        var result = ""
        var first = true

        for ((i, c) in coefficients.withIndex()) {
            val sign = c.sign
            if (sign != 0.0) {
                val number = c.absoluteValue.let {
                    when {
                        i == coefficients.lastIndex -> it.toString()
                        it != 1.0 -> it.toString()
                        else -> ""
                    }
                }
                val x = when (val v = n - i) {
                    0 -> ""
                    1 -> "x"
                    else -> "x" + v.toString().toPow()
                }
                val minus = sign == -1.0

                result += when {
                    first && minus ->  "-"
                    first -> ""
                    minus -> " - "
                    else -> " + "
                } + number + x

                first = false
            }
        }

        return result.takeIf { it.isNotEmpty() } ?: "0"
    }

    fun toFullString(): String {
        var result = ""

        for ((i, c) in coefficients.withIndex()) {
            result += when {
                i == 0 && c < 0 -> "-"
                i == 0 -> ""
                c < 0 -> " - "
                else -> " + "
            } + c.absoluteValue.toString() + ("x" + (n - i).toString().toPow())
        }

        return result
    }

    private fun String.toPow(): String {
        return this
            .replace('0', '⁰')
            .replace('1', '¹')
            .replace('2', '²')
            .replace('3', '³')
            .replace('4', '⁴')
            .replace('5', '⁵')
            .replace('6', '⁶')
            .replace('7', '⁷')
            .replace('8', '⁸')
            .replace('9', '⁹')
    }

    private operator fun get(i: Int) = coefficients[i]

    operator fun unaryMinus() = Polynomial(coefficients.map { -it })

    operator fun unaryPlus() = this

    operator fun plus(a: Double) = Polynomial(coefficients.take(n) + (this[n] + a))

    operator fun minus(a: Double) = plus(-a)

    operator fun times(a: Double) = Polynomial(coefficients.map { it * a })

    operator fun div(a: Double) = Polynomial(coefficients.map { it / a })

    operator fun plus(that: Polynomial): Polynomial {
        val n = max(this.n, that.n) + 1
        val coefficients = MutableList(n) { 0.0 }

        for (i in this.coefficients.indices) {
            coefficients[n - this.n - 1 + i] = this[i]
        }

        for (i in that.coefficients.indices) {
            coefficients[n - that.n - 1 + i] += that[i]
        }

        return Polynomial(coefficients)
    }

    operator fun minus(that: Polynomial): Polynomial {
        val n = max(this.n, that.n) + 1
        val coefficients = MutableList(n) { 0.0 }

        for (i in this.coefficients.indices) {
            coefficients[n - this.n - 1 + i] = this[i]
        }

        for (i in that.coefficients.indices) {
            coefficients[n - that.n - 1 + i] -= that[i]
        }

        return Polynomial(coefficients)
    }

    operator fun times(that: Polynomial): Polynomial {
        val n = this.n + that.n + 1
        val coefficients = MutableList(n) { 0.0 }

        for (i in 0..this.n)
            for (j in 0..that.n)
                coefficients[i + j] += this[i] * that[j]

        return Polynomial(coefficients)
    }

//    data class LongDivisionResult(val quotient: Polynomial, val remainder: Polynomial)
//
//    operator fun div(that: Polynomial): LongDivisionResult {
//        var q = Polynomial(0.0)
//        var r = Polynomial(this.coefficients)
//        val zeroP = Polynomial(0.0)
//
//        while (zeroP != r && r.n >= that.n) {
//            val rc = r.coefficients.find { it != 0.0 } ?: 0.0
//            val dc = that.coefficients.find { it != 0.0 } ?: 0.0
//
//            val tc = rc / dc
//            val te = r.n - that.n
//            val t = Polynomial(tc, *DoubleArray(te))
//
//            q += t
//            r -= t * that
//        }
//
//        return LongDivisionResult(q, r)
//    }

    fun pow(e: Int) = when {
        e == 0 -> One
        n == -1 -> Zero
        e == 1 -> this
        else -> {
            var poly = this
            repeat(e - 1) { poly *= this }
            poly
        }
    }

    operator fun component1() = coefficients.component1()
    operator fun component2() = coefficients.component2()
    operator fun component3() = coefficients.component3()
    operator fun component4() = coefficients.component4()
    operator fun component5() = coefficients.component5()

    override fun equals(other: Any?) =
        if (other !is Polynomial) false else simplify().coefficients == other.simplify().coefficients

    override fun hashCode() = simplify().coefficients.hashCode()

    fun toBezierControls(): List<Double> {
        var list = coefficients.asReversed().mapIndexed { i, c -> c / binomial(n, i) }
        val controls = mutableListOf<Double>()
        while (list.size != 1) {
            controls += list.first()
            list = list.zipWithNext { a, b -> a + b }
        }
        controls += list.first()

        return controls
    }
}

fun List<Polynomial>.sum() = reduceOrNull { a, b -> a + b } ?: Zero

fun List<Polynomial>.product() = reduceOrNull { a, b -> a * b } ?: One


fun solve(coefficients: List<Double>): List<Double> {
    if (coefficients.firstOrNull() == 0.0) {
        return solve(coefficients.dropWhile { it == 0.0 })
    }

    if (coefficients.size <= 1) {
        return listOf()
    }

    if (coefficients.size == 2) {
        val a = coefficients[0]
        val b = coefficients[1]
        return if (a != 0.0) listOf(-b / a) else listOf()
    }

    if (coefficients.size == 3) {
        val a = coefficients[0]
        val b = coefficients[1] / a
        val c = coefficients[2] / a

        val discrim = b * b - 4 * c
        return when {
            discrim approx 0.0 -> {
                val r = -b / 2
                listOf(r, r)
            }
            discrim > 0 -> {
                val e = sqrt(discrim)
                listOf((-b + e) / 2, (-b - e) / 2)
            }
            else -> listOf()
        }
    }

    /*if (coefficients.size == 4) {
        val a = coefficients[0]
        val b = coefficients[1] / a
        val c = coefficients[2] / a
        val d = coefficients[3] / a
        val sqrt3 = sqrt(3.0)

        val f = (3 * c - b * b) / 3
        val g = (2 * b * b * b - 9 * b * c + 27 * d) / 27
        val h = g * g / 4 + f * f * f / 27

        return when {
            f approx 0.0 && g approx 0.0 && h approx 0.0 -> {
                // all roots are real and equal
                val root = cbrt(d)
                listOf(
                    root,
                    root,
                    root
                )
            }

            h <= 0 -> {
                // all 3 roots are real
                val i = sqrt(g * g / 4 - h)
                val j = cbrt(i)
                val k = acos(-(g / (2 * i)))
                val l = j * -1
                val m = cos(k / 3)
                val n = sqrt3 * sin(k/3)
                val p = -(b / 3)

                listOf(
                    2 * j * m - b / 3,
                    l * (m+n) + p,
                    l * (m-n) + p
                )
            }

            else -> {
                // one real and two imaginary roots
                val r = -g / 2 + sqrt(h)
                val s = cbrt(r)
                val t = -g / 2- sqrt(h)
                val u = cbrt(t)

                listOf((s + u) - (b / 3))
            }
        }
    }*/

    /*if (coefficients.size == 5) {
        val a = coefficients[0]
        val b = coefficients[1] / a
        val c = coefficients[2] / a
        val d = coefficients[3] / a
        val e = coefficients[4] / a

        val z0 = b / 4
        val b2 = b * b
        val p = -3 * b2 / 8 + c
        val q = b * b2 / 8 - 1.0 / 2 * b * c + d
        val r = -3.0 / 256 * b2 * b2 + c * b2 / 16 - b * d / 4 + e

        val y = solve(listOf(8.0, -4 * p, -8 * r, 4 * r * p - q * q)).last()

        val a0 = sqrt(-p + 2 * y)
        val b0 = if (a0 == 0.0) y * y - r else -q / 2 / a0

        val r1r2 = solve(listOf(1.0, a0, y + b0))
        val r3r4 = solve(listOf(1.0, -a0, y - b0))
        return (r1r2 + r3r4).map { it - z0 }
    }*/

    return JenkinsTraub.findRoots(coefficients)
}


internal fun det(matrix: List<List<Polynomial>>): Polynomial {
    val n = matrix.size
    if (n == 1) {
        return matrix[0][0]
    }

    if (n == 2) {
        return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]
    }

    if (n == 3) {
        val x = matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]
        val y = matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0]
        val z = matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]
        return x * matrix[0][0] - y * matrix[0][1] + z * matrix[0][2]
    }

//    val sum1 = List(n) { i -> List(n) { j -> matrix[j        ][(i + j) % n] }.product() }.sum()
//    val sum2 = List(n) { i -> List(n) { j -> matrix[n - j - 1][(i + j) % n] }.product() }.sum()
//    val ret = sum1 - sum2

    var det = Zero
    val subN = n - 1
    val temp = List(subN) { MutableList(subN) { Zero } }
    for (col in 0 until n) {
        for (row in 1 until n) {
            for (subCol in 0 until col) {
                temp[row - 1][subCol] = matrix[row][subCol]
            }
            for (subCol in col + 1 until n) {
                temp[row - 1][subCol - 1] = matrix[row][subCol]
            }
        }

        det += matrix[0][col] * det(temp) * (-1.0).pow(col)
    }

    return det
}